This notebook will demonstrate how to use Sentinel-1 backscatter to detect the presence of water.
Sentinel-1 carries a synthetic aperture radar (SAR) sensor, which can observe our planet even through cloud cover. Sentinel-1's SAR sensor is an active sensor; it emits light in the microwave range of the electormagnetic spectrum, and measures the amount that returns, known as backscatter. Smooth surfaces, like calm water, have very low backscatter - most of the light bounces away from the sensor (known as specular reflection). For rough surfaces, like dirt or vegetation, light is scattered in all directions, producing small backscatter signals (known as diffuse backscatter). Large, human-made structures, which feature both vertical and horizonal smooth surfaces, produce large backscatter signals (known as double bounce). As such, the intensity of Sentinel-1 backscatter can distinguish water from land, as shown in the image below:
In the image, the river appears dark (specular reflection), with the urban area appearing very bright (double bounce). For more information, see the article from GIS Geography on SAR.
As of this workshop, Sentinel-1 backscatter is not yet available over South East Asia. Instead, this notebook will demonstrate how to load Sentinel-1 backscatter data over Africa, using Digital Earth Africa's STAC metadata. The Open Data Cube provides odc-stac, a useful library for loading data from STAC into an xarray, helping you work with it in the EASI platform. We will use this library in the workshop.
Important: If you want to work with Sentinel-1 as part of your project, please let us know and we can explore options to make data available to you.
The notebook contains the following steps:
%matplotlib inline
# import required packages
import os
import matplotlib.pyplot as plt
import numpy as np
from dask.diagnostics import ProgressBar
from datacube.utils.cog import write_cog
from dea_tools.plotting import rgb, display_map
from odc.stac import configure_rio, stac_load
from pystac_client import Client
from scipy.ndimage import uniform_filter, variance
from skimage.filters import threshold_minimum
# make a results directory to store output
if not os.path.exists("results"):
os.makedirs("results")
Digital Earth Africa data is stored on S3 in Cape Town, Africa. To load the data, we must configure rasterio with the appropriate AWS S3 endpoint. This is done with the odc.stac.configure_rio function.
The configuration below must be used when loading any Digital Earth Africa data through the STAC API.
# set AWS configuration with specific endpoint
configure_rio(
cloud_defaults=True,
aws={"aws_unsigned": True},
AWS_S3_ENDPOINT="s3.af-south-1.amazonaws.com",
)
Be patient, as this will take a few minutes
# Open the stac catalogue
catalog = Client.open("https://explorer.digitalearth.africa/stac")
catalog
| ID: DEAfrica_data |
| Title: Digital Earth Africa |
| Description: Configure stac endpoint information in your Explorer `settings.env.py` file |
| type: Catalog |
| conformsTo: ['https://api.stacspec.org/v1.0.0-beta.2/core', 'https://api.stacspec.org/v1.0.0-beta.2/item-search'] |
| ID: Arrivals |
| Title: Dataset Arrivals |
| Description: The most recently added Items to this index |
| title: Dataset Arrivals |
| type: Collection |
| properties: {} |
| providers: [] |
| stac_extensions: [] |
| ID: 2e7917da-8db5-5380-84dc-5716ca447afa |
| Bounding Box: [8.329440145578067, 37.842671646059465, 10.948502670082485, 39.95355032947056] |
| Datetime: 2022-09-29 10:00:59.906526+00:00 |
| title: fc_ls_192033_2022-09-29 |
| gsd: 30.0 |
| created: 2022-10-05T03:48:52.037663Z |
| proj:epsg: 32632 |
| platform: landsat-8 |
| landsat:rmse: 7.097 |
| odc:producer: digitalearthafrica.org |
| instruments: ['oli', 'tirs'] |
| eo:cloud_cover: 39.58 |
| view:sun_azimuth: 154.30324656 |
| landsat:rmse_x: 4.804 |
| landsat:rmse_y: 5.225 |
| landsat:wrs_row: 33 |
| odc:file_format: GeoTIFF |
| odc:region_code: 192033 |
| constellation: Landsat |
| view:sun_elevation: 45.51499374 |
| landsat:scene_id: LC81920332022272LGN00 |
| landsat:wrs_path: 192 |
| odc:product_family: fc |
| odc:dataset_version: 1.1.0 |
| odc:naming_conventions: deafrica |
| landsat:collection_number: 2 |
| landsat:collection_category: T1 |
| proj:shape: [7821, 7691] |
| proj:transform: [30.0, 0.0, 439485.0, 0.0, -30.0, 4423215.0, 0.0, 0.0, 1.0] |
| datetime: 2022-09-29T10:00:59.906526Z |
| cubedash:region_code: 192033 |
| stac_extensions: ['https://stac-extensions.github.io/eo/v1.0.0/schema.json', 'https://stac-extensions.github.io/projection/v1.0.0/schema.json', 'https://stac-extensions.github.io/view/v1.0.0/schema.json'] |
| https://stac-extensions.github.io/eo/v1.0.0/schema.json |
| https://stac-extensions.github.io/projection/v1.0.0/schema.json |
| https://stac-extensions.github.io/view/v1.0.0/schema.json |
| href: s3://deafrica-services/fc_ls/1-1-0/192/033/2022/09/29/fc_ls_192033_2022-09-29_bs.tif |
| Title: bs |
| Media type: image/tiff; application=geotiff; profile=cloud-optimized |
| Roles: ['data'] |
| Owner: |
| eo:bands: [{'name': 'bs'}] |
| proj:epsg: 32632 |
| proj:shape: [7821, 7691] |
| proj:transform: [30.0, 0.0, 439485.0, 0.0, -30.0, 4423215.0, 0.0, 0.0, 1.0] |
| href: s3://deafrica-services/fc_ls/1-1-0/192/033/2022/09/29/fc_ls_192033_2022-09-29_pv.tif |
| Title: pv |
| Media type: image/tiff; application=geotiff; profile=cloud-optimized |
| Roles: ['data'] |
| Owner: |
| eo:bands: [{'name': 'pv'}] |
| proj:epsg: 32632 |
| proj:shape: [7821, 7691] |
| proj:transform: [30.0, 0.0, 439485.0, 0.0, -30.0, 4423215.0, 0.0, 0.0, 1.0] |
| href: s3://deafrica-services/fc_ls/1-1-0/192/033/2022/09/29/fc_ls_192033_2022-09-29_ue.tif |
| Title: ue |
| Media type: image/tiff; application=geotiff; profile=cloud-optimized |
| Roles: ['data'] |
| Owner: |
| eo:bands: [{'name': 'ue'}] |
| proj:epsg: 32632 |
| proj:shape: [7821, 7691] |
| proj:transform: [30.0, 0.0, 439485.0, 0.0, -30.0, 4423215.0, 0.0, 0.0, 1.0] |
| href: s3://deafrica-services/fc_ls/1-1-0/192/033/2022/09/29/fc_ls_192033_2022-09-29_npv.tif |
| Title: npv |
| Media type: image/tiff; application=geotiff; profile=cloud-optimized |
| Roles: ['data'] |
| Owner: |
| eo:bands: [{'name': 'npv'}] |
| proj:epsg: 32632 |
| proj:shape: [7821, 7691] |
| proj:transform: [30.0, 0.0, 439485.0, 0.0, -30.0, 4423215.0, 0.0, 0.0, 1.0] |
| href: s3://deafrica-services/fc_ls/1-1-0/192/033/2022/09/29/fc_ls_192033_2022-09-29_thumbnail.jpg |
| Title: Thumbnail image |
| Media type: image/jpeg |
| Roles: ['thumbnail'] |
| Owner: |
| href: s3://deafrica-services/fc_ls/1-1-0/192/033/2022/09/29/fc_ls_192033_2022-09-29.sha1 |
| Media type: text/plain |
| Roles: ['metadata'] |
| Owner: |
| href: s3://deafrica-services/fc_ls/1-1-0/192/033/2022/09/29/fc_ls_192033_2022-09-29.proc-info.yaml |
| Media type: text/yaml |
| Roles: ['metadata'] |
| Owner: |
| Rel: self |
| Target: https://explorer.digitalearth.africa/stac/collections/fc_ls/items/2e7917da-8db5-5380-84dc-5716ca447afa |
| Media Type: application/json |
ODC Dataset YAML
| Rel: odc_yaml |
| Target: https://explorer.digitalearth.africa/dataset/2e7917da-8db5-5380-84dc-5716ca447afa.odc-metadata.yaml |
| Media Type: text/yaml |
| Rel: collection |
| Target: https://explorer.digitalearth.africa/stac/collections/fc_ls |
ODC Product Overview
| Rel: product_overview |
| Target: https://explorer.digitalearth.africa/product/fc_ls |
| Media Type: text/html |
ODC Dataset Overview
| Rel: alternative |
| Target: https://explorer.digitalearth.africa/dataset/2e7917da-8db5-5380-84dc-5716ca447afa |
| Media Type: text/html |
| Rel: items |
| Target: https://explorer.digitalearth.africa/stac/arrivals/items |
Digital Earth Africa
| Rel: root |
| Target: |
| Media Type: application/json |
| Rel: self |
| Target: https://explorer.digitalearth.africa/stac/arrivals |
| Media Type: application/json |
Digital Earth Africa
| Rel: parent |
| Target: |
| Media Type: application/json |
| ID: 2ba5be73-e8e7-5eed-84b9-b658e347e8e8 |
| Bounding Box: [-20.0, -40.000001192092896, 55.00000111758709, 40.0] |
| Datetime: 1981-01-01 00:00:00+00:00 |
| title: chirps-v2.0_1981.01.01.stac-item |
| proj:epsg: 4326 |
| proj:shape: [1600, 1500] |
| odc:product: rainfall_chirps_daily |
| proj:transform: [0.05000000074505806, 0.0, -20.0, 0.0, -0.05000000074505806, 40.0, 0.0, 0.0, 1.0] |
| odc:file_format: GeoTIFF |
| end_datetime: 1981-01-01T23:59:59Z |
| start_datetime: 1981-01-01T00:00:00Z |
| created: 2021-12-16T05:00:16.622284Z |
| datetime: 1981-01-01T00:00:00Z |
| cubedash:region_code: None |
| stac_extensions: ['https://stac-extensions.github.io/eo/v1.0.0/schema.json', 'https://stac-extensions.github.io/projection/v1.0.0/schema.json'] |
| https://stac-extensions.github.io/eo/v1.0.0/schema.json |
| https://stac-extensions.github.io/projection/v1.0.0/schema.json |
| href: s3://deafrica-input-datasets/rainfall_chirps_daily/1981/01/chirps-v2.0_1981.01.01.tif |
| Title: rainfall |
| Media type: image/tiff; application=geotiff; profile=cloud-optimized |
| Roles: ['data'] |
| Owner: |
| eo:bands: [{'name': 'rainfall'}] |
| proj:epsg: 4326 |
| proj:shape: [1600, 1500] |
| proj:transform: [0.05000000074505806, 0.0, -20.0, 0.0, -0.05000000074505806, 40.0, 0.0, 0.0, 1.0] |
| Rel: self |
| Target: https://explorer.digitalearth.africa/stac/collections/rainfall_chirps_daily/items/2ba5be73-e8e7-5eed-84b9-b658e347e8e8 |
| Media Type: application/json |
ODC Dataset YAML
| Rel: odc_yaml |
| Target: https://explorer.digitalearth.africa/dataset/2ba5be73-e8e7-5eed-84b9-b658e347e8e8.odc-metadata.yaml |
| Media Type: text/yaml |
| Rel: collection |
| Target: https://explorer.digitalearth.africa/stac/collections/rainfall_chirps_daily |
ODC Product Overview
| Rel: product_overview |
| Target: https://explorer.digitalearth.africa/product/rainfall_chirps_daily |
| Media Type: text/html |
ODC Dataset Overview
| Rel: alternative |
| Target: https://explorer.digitalearth.africa/dataset/2ba5be73-e8e7-5eed-84b9-b658e347e8e8 |
| Media Type: text/html |
Digital Earth Africa
| Rel: root |
| Target: |
| Media Type: application/json |
| Rel: self |
| Target: https://explorer.digitalearth.africa/stac |
| Media Type: application/json |
Collections
| Rel: children |
| Target: https://explorer.digitalearth.africa/stac/collections |
| Media Type: application/json |
| description: All product collections |
Arrivals
| Rel: child |
| Target: |
| Media Type: application/json |
| description: Most recently added items |
Item Search
| Rel: search |
| Target: https://explorer.digitalearth.africa/stac/search |
| Media Type: application/json |
alos_palsar_mosaic
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/alos_palsar_mosaic |
| description: ALOS/PALSAR and ALOS-2/PALSAR-2 annual mosaic tiles generated for use in the Data Cube - 25m pixel spacing, WGS84. These tiles are derived from the orignal JAXA mosaics with conversion to GeoTIFF. |
cci_landcover
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/cci_landcover |
| description: ESA Climate Change Initiative Land Cover |
cgls_landcover
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/cgls_landcover |
| description: Copernicus Global Land Service, Land Use/Land Cover at 100 m |
crop_mask
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/crop_mask |
| description: Annual cropland extent map produced by Digital Earth Africa. |
crop_mask_central
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/crop_mask_central |
| description: Annual cropland extent map for Central Africa produced by Digital Earth Africa. |
crop_mask_eastern
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/crop_mask_eastern |
| description: Annual cropland extent map for Eastern Africa produced by Digital Earth Africa. |
crop_mask_indian_ocean
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/crop_mask_indian_ocean |
| description: Annual cropland extent map for Indian Ocean Africa produced by Digital Earth Africa. |
crop_mask_northern
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/crop_mask_northern |
| description: Annual cropland extent map for Northern Africa produced by Digital Earth Africa. |
crop_mask_sahel
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/crop_mask_sahel |
| description: Annual cropland extent map for Sahel Africa produced by Digital Earth Africa. |
crop_mask_southeast
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/crop_mask_southeast |
| description: Annual cropland extent map for Southeast Africa produced by Digital Earth Africa. |
crop_mask_southern
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/crop_mask_southern |
| description: Annual cropland extent map for Southern Africa produced by Digital Earth Africa. |
crop_mask_western
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/crop_mask_western |
| description: Annual cropland extent map for Western Africa produced by Digital Earth Africa. |
dem_cop_30
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/dem_cop_30 |
| description: Copernicus DEM 30 m |
dem_cop_90
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/dem_cop_90 |
| description: Copernicus DEM 90 m |
dem_srtm
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/dem_srtm |
| description: 1 second elevation model |
dem_srtm_deriv
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/dem_srtm_deriv |
| description: 1 second elevation model derivatives |
esa_worldcover
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/esa_worldcover |
| description: ESA World Cover, global 10 m land use/land cover data from 2020. |
fc_ls
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/fc_ls |
| description: Landsat Fractional Cover Observations from Space |
fc_ls_summary_annual
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/fc_ls_summary_annual |
| description: DE Africa Landsat Fractional Cover Percentiles |
gm_ls5_ls7_annual
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/gm_ls5_ls7_annual |
| description: Surface Reflectance Annual Geometric Median and Median Absolute Deviations, Landsat 5 and Landsat 7 |
gm_ls5_ls7_annual_lowres
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/gm_ls5_ls7_annual_lowres |
| description: Surface Reflectance Annual Geometric Median and Median Absolute Deviations, Landsat 5 and Landsat 7. Low resolution version used for visualisations. |
gm_ls8_annual
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/gm_ls8_annual |
| description: Surface Reflectance Annual Geometric Median and Median Absolute Deviations, Landsat 8 |
gm_ls8_annual_lowres
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/gm_ls8_annual_lowres |
| description: Surface Reflectance Annual Geometric Median and Median Absolute Deviations, Landsat 8. Low resolution version used for visualisations. |
gm_ls8_ls9_annual
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/gm_ls8_ls9_annual |
| description: Surface Reflectance Annual Geometric Median and Median Absolute Deviations, Landsat 8 and Landsat 9 |
gm_ls8_ls9_annual_lowres
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/gm_ls8_ls9_annual_lowres |
| description: Surface Reflectance Annual Geometric Median and Median Absolute Deviations, Landsat 8 and Landsat 9. Low resolution version used for visualisations. |
gm_s2_annual
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/gm_s2_annual |
| description: Surface Reflectance Annual Geometric Median and Median Absolute Deviations, Sentinel-2 |
gm_s2_annual_lowres
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/gm_s2_annual_lowres |
| description: Annual Geometric Median, Sentinel-2 - Low Resolution mosaic |
gm_s2_semiannual
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/gm_s2_semiannual |
| description: Surface Reflectance Semiannual Geometric Median and Median Absolute Deviations, Sentinel-2 |
gm_s2_semiannual_lowres
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/gm_s2_semiannual_lowres |
| description: Surface Reflectance Semiannual Geometric Median and Median Absolute Deviations, Sentinel-2. Low resolution version used for visualisations. |
gmw
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/gmw |
| description: Global Mangrove Watch data sourced from the UN Environment Program at https://data.unep-wcmc.org/datasets/45 |
io_lulc
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/io_lulc |
| description: Impact Observatory (ESRI) Landcover Classification 9 Class |
jers_sar_mosaic
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/jers_sar_mosaic |
| description: JERS-1 SAR annual mosaic tiles generated for use in the Data Cube 25m pixel spacing, WGS84. These tiles are derived from the orignal JAXA mosaics with conversion to GeoTIFF. |
ls5_sr
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/ls5_sr |
| description: USGS Landsat 5 Collection 2 Level-2 Surface Reflectance |
ls5_st
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/ls5_st |
| description: USGS Landsat 5 Collection 2 Level-2 Surface Temperature |
ls7_sr
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/ls7_sr |
| description: USGS Landsat 7 Collection 2 Level-2 Surface Reflectance |
ls7_st
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/ls7_st |
| description: USGS Landsat 7 Collection 2 Level-2 Surface Temperature |
ls8_sr
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/ls8_sr |
| description: USGS Landsat 8 Collection 2 Level-2 Surface Reflectance |
ls8_st
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/ls8_st |
| description: USGS Landsat 8 Collection 2 Level-2 Surface Temperature |
ls9_sr
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/ls9_sr |
| description: USGS Landsat 9 Collection 2 Level-2 Surface Reflectance |
ls9_st
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/ls9_st |
| description: USGS Landsat 9 Collection 2 Level-2 Surface Temperature |
nasadem
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/nasadem |
| description: NASADEM from Microsoft's Planetary Computer |
ndvi_anomaly
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/ndvi_anomaly |
| description: Monthly NDVI Anomalies produced by Digital Earth Africa. |
ndvi_climatology_ls
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/ndvi_climatology_ls |
| description: Monthly NDVI Climatologies produced by Digital Earth Africa. |
pc_s2_annual
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/pc_s2_annual |
| description: Surface Reflectance Annual Clear Pixel Count, Sentinel-2 |
rainfall_chirps_daily
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/rainfall_chirps_daily |
| description: Rainfall Estimates from Rain Gauge and Satellite Observations |
rainfall_chirps_monthly
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/rainfall_chirps_monthly |
| description: Rainfall Estimates from Rain Gauge and Satellite Observations |
s1_rtc
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/s1_rtc |
| description: Sentinel 1 Gamma0 normalised radar backscatter |
s2_l2a
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/s2_l2a |
| description: Sentinel-2a and Sentinel-2b imagery, processed to Level 2A (Surface Reflectance) and converted to Cloud Optimized GeoTIFFs |
wofs_ls
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/wofs_ls |
| description: Historic Flood Mapping Water Observations from Space |
wofs_ls_summary_alltime
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/wofs_ls_summary_alltime |
| description: Water Observations from Space Alltime Statistics |
wofs_ls_summary_annual
| Rel: child |
| Target: https://explorer.digitalearth.africa/stac/collections/wofs_ls_summary_annual |
| description: Water Observations from Space Annual Statistics |
Digital Earth Africa
| Rel: root |
| Target: |
| Media Type: application/json |
The specifics of Digital Earth Africa's Sentinel-1 backscatter product are available in the Digital Earth Africa Explorer
The data has five measurements, we will look at two today:
vv is the amount of vertically polarised backscatter received by the sensorvh is the amount of horizontally polarised backscatter received by the sensorPolarisation is the orientation of the light relative to the sensor, and different scattering surfaces can effect the backscatter polarisation. In the case of Sentinel-1, the sensor sends out vertically polarised light (v) and measures the amount of vertically polarised light (vv) or horizontally polarised light (vh) returned.
Digital Earth Africa's Sentinel-1 data (s1_rtc) needs to be configured with additional information about the data type and no-data values. This is provided in the dictionary below:
The configuration dictionary is determined from the product’s definition, available at https://explorer.digitalearth.africa/products/s1_rtc
# VV, VH, AREA = NaN, float32
# angle = 255, uint8
# mask = 0, uint8
# set the collection configuration, setting the products definition
config = {
"s1_rtc": {
"assets": {
"*": {
"data_type": "float32",
"nodata": float("nan"),
"unit": "1",
},
"angle": {
"data_type": "uint8",
"nodata": 255,
"unit": "1",
},
"mask": {
"data_type": "uint8",
"nodata": 0,
"unit": "1",
},
},
}
}
For this example, we'll use a bounding box, but you could come back an update this to use a polygon by following training session 2
When working with STAC, we use the word "collection" to refer to a dataset, rather than "product" when working with the ODC. For more information on the differences, see the odc-stac documentation
# Set a bounding box using latitude and longitude coordinates
# [xmin, ymin, xmax, ymax]
bbox = [-16.34, 12.5699, -16.24, 12.67]
# timeframe
single_date = "2020"
# Set the STAC collections
collections = ["s1_rtc"]
# View area of interest
display_map(x=(bbox[0], bbox[2]), y=(bbox[1],bbox[3]))
The connection to Digital Earth Africa's STAC endpoint provides a search function which returns dataset ids matching the parameters of the search. Data can then be loaded by accessing the corresponding files for each id.
# Build a query with the set parameters
query = catalog.search(bbox=bbox, collections=collections, datetime=single_date)
# Search the STAC catalog for all items matching the query
items = list(query.get_items())
print(f"Found: {len(items):d} datasets")
Found: 31 datasets
In this section we will use the stac_load function to only load data within the defined parameters set above.
This data will be lazy-loaded with dask, meaning the it will not be loaded into memory until required, such as when it is displayed.
crs = "EPSG:6933"
resolution = 20
ds = stac_load(
items,
bands=[
"vv",
"vh",
],
crs=crs,
resolution=resolution,
chunks={},
groupby="solar_day",
stac_cfg=config,
bbox=bbox,
)
# View the Xarray Dataset - loads as dask array
ds
<xarray.Dataset>
Dimensions: (y: 624, x: 484, time: 31)
Coordinates:
* y (y) float64 1.604e+06 1.604e+06 ... 1.591e+06 1.591e+06
* x (x) float64 -1.577e+06 -1.577e+06 ... -1.567e+06 -1.567e+06
spatial_ref int32 6933
* time (time) datetime64[ns] 2020-01-04T19:17:34.083167 ... 2020-12...
Data variables:
vv (time, y, x) float32 dask.array<chunksize=(1, 624, 484), meta=np.ndarray>
vh (time, y, x) float32 dask.array<chunksize=(1, 624, 484), meta=np.ndarray># Load the data in memory using dask's .compute() function
with ProgressBar():
ds = ds.compute()
[########################################] | 100% Completed | 11.15 s
The vv and vh bands can be viewed individually, or as a RGB style image through combination with their ratio as the third band.
# Selecting a few images from the loaded S1 to visualise
timesteps = [2, 4, 6, 9, 11]
# Plot VV polarisation for specific timeframe
ds.vv.isel(time=timesteps).plot(cmap="Greys_r", robust=True, col="time", col_wrap=5);
# Plot VH polarisation for specific timeframe
ds.vh.isel(time=timesteps).plot(cmap="Greys_r", robust=True, col="time", col_wrap=5);
For the RGB visualization below, the ratio between VH and VV is added as a third measurement band.
# VH/VV is a potentially useful third feature after VV and VH
ds["vh/vv"] = ds.vh / ds.vv
# median values are used to scale the measurements so they have a similar range for visualization
med_s1 = ds.median()
# plotting an RGB image for selected timesteps
rgb(
ds[["vv", "vh", "vh/vv"]] / med_s1,
bands=["vv", "vh", "vh/vv"],
index=timesteps,
col_wrap=5,
);
The next three cells will export the first observation (time=0) to a GeoTiff.
If you want to export multiple observations, refer to training session 2.
da = ds.to_array()
rgb_ga = da.isel(time=0)
write_cog(geo_im=rgb_ga, fname="results/sar_rgb.tif", overwrite=True)
PosixPath('results/sar_rgb.tif')
A common processing step when working with SAR data is to remove , which causes SAR data to have a grainy appearence. Radar observations appear speckly) due to random interference of coherent signals from target scatters. The speckle noise can be reduced by averaging pixel values over an area or over time. However, averaging over a fixed window smoothes out real local spatial variation and leads to reduced spatial resolution. An adaptive approach that takes into account local homogeneity is therefore preferred.
Below, we display the unfiltered data, then apply the Lee filter, one of the popular adaptive speckle filters.
# images appear smoother after speckle filtering
fig, ax = plt.subplots(1, 2, figsize=(15, 5))
ds["vv"].isel(time=3).plot(ax=ax[0], robust=True)
ds["vh"].isel(time=3).plot(ax=ax[1], robust=True)
ax[0].set_title("VV")
ax[1].set_title("VH")
plt.tight_layout();
# defining a function to apply lee filtering on S1 image
def lee_filter(da, size):
"""
Apply lee filter of specified window size.
Adapted from https://stackoverflow.com/questions/39785970/speckle-lee-filter-in-python
"""
img = da.values
img_mean = uniform_filter(img, (size, size))
img_sqr_mean = uniform_filter(img**2, (size, size))
img_variance = img_sqr_mean - img_mean**2
overall_variance = variance(img)
img_weights = img_variance / (img_variance + overall_variance)
img_output = img_mean + img_weights * (img - img_mean)
return img_output
Now that we’ve defined the filter, we can run it on the VV and VH data. You might have noticed that the function takes a size argument. This will change how blurred the image becomes after smoothing. We’ve picked a default value for this analysis, but you can experiement with this if you’re interested.
# The lee filter above doesn't handle null values
# We therefore set null values to 0 before applying the filter
valid = ds.where(~np.isinf(ds))
ds = ds.where(valid, 0)
# Create a new entry in dataset corresponding to filtered VV and VH data
ds["filtered_vv"] = ds.vv.groupby("time").apply(lee_filter, size=7)
ds["filtered_vh"] = ds.vh.groupby("time").apply(lee_filter, size=7)
# Null pixels should remain null
ds["filtered_vv"] = ds.filtered_vv.where(valid.vv, np.nan)
ds["filtered_vh"] = ds.filtered_vh.where(valid.vh, np.nan)
# images appear smoother after speckle filtering
fig, ax = plt.subplots(1, 2, figsize=(15, 5))
ds["filtered_vv"].isel(time=3).plot(ax=ax[0], robust=True)
ds["filtered_vh"].isel(time=3).plot(ax=ax[1], robust=True)
ax[0].set_title("VV")
ax[1].set_title("VH")
plt.tight_layout();
Similar to training session 2, it is possible to view the histograms for the vv and vh bands and identify an automatic threshold that seperates two key features in the landscape: water, which has very low backscatter, and land, which has higher backscatter.
When working with radar backscatter, it is common to work with the data in units of decibels (dB), rather than digital number (DN) as measured by the satellite. To convert from DN to dB, we use the following formula:
$$\text{dB} = 10 \times \log_{10}(\text{DN})$$# convert digital numbers to dB
ds["filtered_vv"] = 10 * np.log10(ds.filtered_vv)
ds["filtered_vh"] = 10 * np.log10(ds.filtered_vh)
# histogram analysis for S1
fig = plt.figure(figsize=(12, 3))
ds.filtered_vh.plot.hist(bins=1000, label="VH filtered")
ds.filtered_vv.plot.hist(bins=1000, label="VV filtered", alpha=0.5)
plt.xlim(-40, -1)
plt.legend()
plt.xlabel("DN values in(dB)")
plt.title("Comparison of Lee filtered VH and VV polarisation bands");
The vv band has a clear separation between the two classes. There is a smaller number of low value pixels, peaking around -30 dB (likely water), and a larger number of high value pixels, peaking around -15 dB (likely land).
The histogram for VH backscatter shows a bimodal distribution with low values over water and high values over land. The VV histogram has multiple peaks and less obvious seperation between water and land.
We therefore build a classifier based on VH backscatter. We choose a threshold to separate land and water: pixels with values below the threshold are water, and pixels with values above the threshold are not water (land).
There are several ways to determine the threshold. Here, we use the threshod_minimum function implemented in the skimage package to determine the threshold from the VH histogram automatically. This method computes the histogram for all backscatter values, smooths it until there are only two maxima and find the minimum in between as the threshold.
# Convert to numpy ndarray, remove nan values, and calculate automatic threshold
vh_new = ds["filtered_vh"].values
vh_new = vh_new[~np.isnan(vh_new)]
threshold_vh = threshold_minimum(vh_new)
print(threshold_vh)
-25.661585
# visualise the threshold
fig, ax = plt.subplots(figsize=(15, 3))
ds.filtered_vh.plot.hist(bins=1000, label="VH filtered")
plt.xlim(-40, -5)
ax.axvspan(xmin=-40.0, xmax=threshold_vh, alpha=0.25, color="green", label="Water")
ax.axvspan(xmin=threshold_vh, xmax=-5, alpha=0.25, color="red", label="Not Water")
plt.legend()
plt.xlabel("VH (dB)")
plt.title("Effect of the classifier")
plt.show()
# define the classifier, values lower than the threshold are water pixels, return dataset containing water pixels
def S1_water_classifier(da, threshold=threshold_vh):
water_data_array = da < threshold
return water_data_array.to_dataset(name="s1_water")
# return classified data product
ds["water"] = S1_water_classifier(ds.filtered_vh).s1_water
ds
<xarray.Dataset>
Dimensions: (time: 31, y: 624, x: 484)
Coordinates:
* y (y) float64 1.604e+06 1.604e+06 ... 1.591e+06 1.591e+06
* x (x) float64 -1.577e+06 -1.577e+06 ... -1.567e+06 -1.567e+06
spatial_ref int32 6933
* time (time) datetime64[ns] 2020-01-04T19:17:34.083167 ... 2020-12...
Data variables:
vv (time, y, x) float32 0.05976 0.04359 0.02704 ... 0.2923 0.1862
vh (time, y, x) float32 0.005369 0.004237 ... 0.01821 0.01918
vh/vv (time, y, x) float32 0.08984 0.09722 0.1888 ... 0.06229 0.103
filtered_vv (time, y, x) float32 -12.21 -12.32 -12.53 ... -7.404 -7.758
filtered_vh (time, y, x) float32 -22.11 -22.5 -22.74 ... -19.43 -19.19
water (time, y, x) bool False False False False ... False False FalseWe can now view the image with our classification. The classifier returns either True or False for each pixel. To investigate the landscape, we want to check which pixels are always water and which are always land. Conveniently, Python encodes True = 1 and False = 0. If we plot the average classified pixel value, pixels that are always water will have an average value of 1 and pixels that are always land will have an average of 0. Pixels that are sometimes water and sometimes land will have an average between these values.
The following cell plots the average classified pixel value over time. What can you learn about the landscape from this summary?
# Plot the mean of each classified pixel value
plt.figure(figsize=(8, 7))
ds.water.mean(dim="time").plot(cmap="RdBu")
plt.title("Average classified pixel value");
We can also see if the standard deviation of each pixel in time is a reasonable way to learn about the landscape. Similar to how we calculated and plotted the mean above, you can calculate and plot the standard deviation by using the std function in place of the mean function. The standard deviation gives us an idea of how variable a pixel has been over the entire period of time that we looked at.
If you’d like to see the results using a different colour-scheme, you can also try substituting cmap="Greys" or cmap="Blues" in place of cmap="viridis".
plt.figure(figsize=(8, 7))
ds.water.std(dim="time").plot(cmap="viridis")
plt.title("Standard deviation of classified pixel values");
Make a copy of this notebook so that you can keep this one as a reference. Some things you might like to try: